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Towards  Aspect  Invariant  Feature  Sets  for 
Characterizing  Three  Dimensional  Objects 

Abstract 
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The  problem  of  characterizing  three  dimensional  objects  with 


a  minimum  number  of  feature  sets  is  addressed.  Previous  research  has 


been  plagued  by  features  which  are  a  function  of  aspect  angle  and  so 
efforts  have  concentrated  on  characterizing  an  object  with  several 
hundred  prototype  feature  sets.  In  this  work  it  is  demonstrated  that 
in  many  cases,  the  silhouette  obtained  from  one  view  of  the  object  can 
be  derived  from  a  linear  transformation  of  the  silhouette  from  another 


view.  As  a  result  of  this  relationship  a  single  set  of  moments  which 
is  invariant  to  such  a  general  linear  transformation  can  be  used  to 


characterize  many  views  of  the  same  object  and  hence  the  number  of 
prototvpe  feature  sets  required  to  specify  an  object  is  reduced.  In 
addition,  it  is  demonstrated  that  for  some  objects  it  is  advantageous 
to  partition  the  object  into  regions  in  order  to  find  the  region  of 
the  object  which  is  least  dependent  upon  aspect  angle. 
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Fhe  problem  of  characterizing  three  dimensional  objects  with 
a  minimum  number  of  feature  sets  is  addressed.  Previous  research 
has  been  plagued  by  features  which  are  a  function  of  aspect  angle 
and  so  efforts  have  concentrated  on  characterizing  an  object  with 
several  hundred  prototype  feature  sets.  In  this  work  it"  is 
demonstrated  that  in  many  cases,  the  silhouette  obtained  from  one 
view  of  the  object  can  be  derived  from  a  linear  transformation  of 
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INTRODUCTION 

For  quite  some  time,  the  recognition  of  hand  printed 
characters  has  received  much  attention  and  the  methods  developed 
have  achieved  significant  success.  In  recent  years,  considerable 
interest  has  developed  in  the  recognition  of  three  dimensional 
objects.  The  research  has  been  largely  motivated  by  computer  aided 
assembly  and  inspection  for  manufacturing  and  air  traffic  control 
problems.  In  the  first  case,  the  object  is  generally  oriented  so 
that  recognition  or  inspection  can  be  accomplished  in  two  dimensions 
(1,3, 4,6).  Some  attempts  have  been  made  to  use  multiple  views  of 
the  object  (2,5)  but  these  methods  are  not  applicable  to  air  traffic 
control  problems.  In  this  case,  information  in  three  dimensions  is 
necessary  for  recognition  but  only  one  view  is  available.  This  makes 
object  recognition  more  difficult  since  parts  of  the  object  may  be 
obscured  by  other  parts  of  the  object.  Thus,  the  shape  of  the  sil¬ 
houette  seen  by  the  camera  may  be  dependent  upon  the  viewing  or  aspect 
angle.  The  orientation  of  the  object  may  further  alter  the  shape  of 
the  silhouette  through  rotation,  translation,  and  size  change  (simili¬ 
tude)  . 

From  previous  work  in  character  recognition,  it  is  well  known 
that  features,  called  moment  invariants  (9-11),  can  be  used  to  describe 
the  shape  of  a  two  dimensional  image.  They  are  invariant  to  rotation, 
translation,  and  similitude.  These  features  have  been  directly  applied 
to  the  recognition  of  three  dimensional  objects  (12-16).  In  (12), 

Dudani  extracts  moment  invariant  features  from  aircraft  shapes  using 
over  500  different  aspect  angles  for  each  aircraft  type.  For  the 
six  different  type  of  aircraft  used,  the  complete  training  sample  set 
contained  over  3000  live  images.  To  classify  an  unknown  aircraft,  a 
modified  K-nearest  neighbor  classifier  was  used  to  find  the  ten  nearest 
neighbors  among  the  3000  element  training  set.  A  high  degree  of  correct 
recognition  was  reported  and  is  probably  the  most  successful  work  to 
date.  The  major  disadvantages  of  the  method  are  the  computational 
load,  storage  requirements,  and  the  large  number  of  images  required  to 
characterize  an  object.  A  similar  aproach  to  recognizing  ships  was 
taken  by  Smith  and  Wright  (13).  Again,  moments  invariant  to  rotation, 
translation,  and  similitude  were  used  to  determine  ship  type.  Here, 
the  aspect  problem  was  ingored  by  permitting  only  top  views  of  ships 
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to  be  admissible  images.  Another  study  for  recognizing  tanks  (14, 

17)  was  done  and  the  methodology  is  the  same  as  that  of  Dudani. 

A  comprehensive  survey  of  automatic  recognition  of  three  dimen¬ 
sional  objects  using  one  optical  sensor  can  be  found  in  McGhee  (15). 

Here,  possible  features  are  discussed  for  pattern  recognition  classi¬ 
fiers,  such  as  the  slope  code  formulation,  Fourier  descriptors  (18), 
and  \cment  invariants.  In  addition,  some  other,  less  successful, 
techniques  for  recognition  of  three  dimensional  objects  are  presented. 

In  all  of  the  above  methods  for  recognizing  three  dimensional 
objects,  features  are  used  which  are  a  function  of  aspect  angle  and 
so  characterizing  an  object  typically  requires  several  hundred  feature 
sets.  Such  an  approach  does  not  consider  the  fact  that  some  regions  of 
the  object  are  informationally  richer  than  others  (20).  In  this  work, 
a  region  which  is  rich  in  information  is  one  that  has  features  which 
are  aspect  invariant  while  all  other  regions  are  considered  to  be 
ambiguous  and  are  not  useful  for  classifying  the  object.  This  idea 
is  consistent  with  the  conjecture  that  shape  recognition  is  a  hierarchical 
process  (19) • 

In  this  work  an  attempt  is  made  to  characterize  three  dimensional 
objects  with  a  minimum  number  of  feature  sets.  Here  it  is  demonstrated 
that  for  many  cases,  the  silhouette  obtained  from  one  view  of  the  object 
can  be  derived  from  a  linear  transformation  of  the  silhouette  from 
another  view.  As  a  result  of  this  relationship,  a  set  of  moments  which 
is  invariant  to  such  a  general  linear  transformation  can  be  used  to 
characterize  many  views  of  the  same  object  and  hence  the  number  of 
feature  sets  required  to  specify  an  object  is  reduced. 

The  next  section  discusses  a  set  of  features  which  are  invariant 
under  a  general  linear  transformation.  This  is  followed  by  some  experi¬ 
ments  with  some  common  solid  objects  in  which  features  are  computed 
for  various  aspect  angles.  Next,  aircraft  shapes  are  considered  and 

features  are  extracted  for  various  regions  of  the  aircraft  at  numerous 
aspect  angles.  Finally,  suggestions  for  future  research  are  presented. 
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Moments  Invariant  Under  General  Linear  ^ 

Tranformat Ions 

Let  the  two  dimensional  (p+q)th  order  moments  of  a  discrete 
density  distribution  fCx^y^)  be  defined  as 

,  N 

=  — fj—  (x1-xi)P  (yi-Yi)  fCXi.Yi)  (1) 

P,q  =  0,1,2.  .  . 

where  x  is  the  mean  of  the  coordinate  x^  and  y^  is  the  mean  of  the 
coordinate  y^.  For  a  solid  silhouette  discretized  to  a  matrix  of 
zeros  and  ones  the  distribution  f(x^,y^)  =  1  for  a  point  contained 
within  the  silhouette  and  zero  otherwise.  It  can  thdn  be  shown  that 
moments  of  all  orders  exist  and  that  the  sequence  XM  }  uniquely 
determined  by  fCx^y^;  and  conversely,  f(x1,y1)  is  uniquely  determined 

by  (H)- 

For  the  general  linear  transformation 


u  _  a  b  x 


LVJ  Lc  tl  LyJ 

there  exists  four  absolute  moment  invariants  given  by 

R1  =  ^M20  M02  "  M11  /  M003  (3) 
r2  =  ((m30  mQ3  -  m21  m12)2  -  it  (m3Q  m12  -  m212)  (m21  mQ3  -  m122))/m 

R3  =  ('M20^M21  M03  "  M02  ^  ~  M11^M30  M03  “  M21  M12^  + 

M02(M30  M12  '  M212))/  M005  (5) 

and 

R4  =(M302  M023  ~  6M30  M21  M11  ^  +  6  M30  M12  M02  (2  Mll2  "  M20  M02) 

+  M30  M03  ^6M20  M11  M02  “  8  Mll3^ 

+  9  M21  M20  MQ2  -  18  M21  M12  M20  mq2  +  6M21  MQ3  M20  ( 2MX1 

M20  M02^ 

+  9  M12  M20  MQ2  -  6M12  MQ3  Mi;l  M20  +  MQ3  M203  )/Mqo7  (6) 


Moments  that  are  invariant  to  size  change,  rotation  and  translation 
have  been  used  extensively  for  characterizing  shape  (11-16)  and  they 
are  clearly  a  function  of  aspect  angle  (15).  Moments  that  are  invariant 


I 


5 


to  a  linear  transformation  can  reduce  this  dependency  on  aspect 
angle.  For  example,  consider  a  cube  centered  about  the  origin  in 
an  xyz  coordinate  system  where  <j>  represents  rotation  about  the  x 
axis,  9  about  the  y  axis  and  is  rotation  about  the  z  axis.  Initially, 
4)  =  9  =  ¥  =  0 .  Rotating  the  object  through  any  angle  9  will  not 
change  the  value  of  R^,R2,R2»  or  R ^  since  the  silhouette  obtained 
after  rotation  through  9  given  by  the  coordinates  (i^,v^)  can  be 
obtained  from  the  original  silhouette  by 


a 

0 

v 

i 

0 

1 

h 

On  the  other  hand,  moments  which  are  invariant  to  size  change, 
rotation  (this  only  includes  rotation  about  an  axis  perpendicular 
to  the  original  viewing  plane)  and  translation  will  change  after 
the  rotation  through  9.  Consequently,  moments  invariant  to  a  linear 
transformation  are  able  to  characterize  more  regions  of  the  object. 
The  net  result  is  a  reduction  In  the  number  of  feature  sets  needed 
to  characterize  the  object. 


Investigating  Aspect  Invariance  for  Some  Simple  Objects 

The  features  R^,  R^,  R^  and  R^  were  tested  on  a  sphere,  cube, 
cylinder  and  a  cone.  All  images  were  binary  images  and  were  simulated 
on  an  IBM  370/158. 

Let  the  object  reside  in  a  coordinate  system  xyz  with  the  ^prigin 
at  or  close  to  the  center  of  gravity  of  the  object.  Initially,  let 
the  three  axes  xyz  be  colinear  with  the  axes  of  the  viewing  or  camera 
coordinate  system,  uvw.  The  viewing  coordinate  system  is  then  rotated 
about  the  y  axis  through  the  angle  9  ,  about  the  z  axis  through  the 
angle  f  and  about  the  x  axis  through  <J>,  in  that  order.  The  transfor¬ 
mation  which  relates  the  two  coordinate  systems  is 


—  — 

-  — 

-  - 

u 

cos  9  cos  V  sin  9  sin  4>  sin  9  cos  4>  + 

X 

cos  9  sin  41  cos  4>  cos  9  sin  ¥  sin  4> 

V 

= 

sin  T  cos  y  cos  4>  -cos  f  sin  <j> 

y 

w 

-sin  9  cos  V  cos  9  sin  <j>  +  cos  9  cos  4>  - 

_  sin  9  sin  'V  cos  <f>  sin  9  sin  V  sin  <f> 

z 

(8) 


The  moments  in  equations  ( 3 ) — ( 6 )  were  then  calculated  from  the 
solid  silhouette  in  the  uv  plane  for  various  orientations  of  the  object 
Now,  the  moments  used  here  will  be  invariant  to  a  rotation  through 
an  angle  Y  since  this  is  just  a  special  case  of  a  linear  transformation 
Furthermore,  for  an  image  that  possesses  symmetry  about  a  plane,  the 
range  of  9  and  <}>  necessary  to  cover  all  distinct  views  in  smaller 
than  the  total  0  to  360  degree  range  (12).  It  can  be  shown  that  the 
significant  range  of  values  is 

-90°  <  9  <90°  (9) 

0°  <  <t>  <  90° 

All  objects  were  represented  in  the  uv  plane  by  an  80  x  80  matrix  of 
ones  and  zeros.  In  the  original  xyz  coordinate  system  the  sphere  can 


be  described  by 

x2  +  y2  +  z2  <  100  (10) 

the  cube  can  be  described  by  the  boundaries 

x  =  ±  5,  y  =  ±5,  z  =  ±5  (ID 

The  cylinder  has  boundaries 

x2  +  y2  £  25,  z  =  ±  5  (12) 

and  the  cone  can  be  described  by 

x2  +  y2  £  25,  z  <  5  -  | x |  ,  z  <  5  -  | y |  tl3) 


A  summary  of  the  results  is  shown  in  Table  1  where  Ri  and  a?  are  the 
mean  and  variance  respectively  of  FL  . 

Although  it  is  not  evident  from  this  table,  was  found  to  be 
very  sensitive  to  noise  or  to  the  slightest  distortion  of  the  object 

p 

This  is  the  reason  for  the  large  variance,  .  As  a  result  of  this 
observation,  was  not  computed  for  the  cube  and  cylinder.  In 
addition,  the  features  R2  and  R^  are  nearly  zero  in  all  cases  except 
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the  cone.  was  the  only  moment  found  to  be  effective  in  measuring 

the  dependency  on  aspect  angle. 

The  quantity  k  refers  to  the  number  of  feature  vectors  that 
are  needed  to  characterize  the  object.  This  was  calculated  from 
the  following  algorithm. 

1.  Set  k=l,  Feat(j,l)=0  for  j=l,2,3>4. 

2.  Compute  R.  for  j=l,2,3,^. 

J 

3.  If  |  R^-Feat  ( 1,  i)  |  +  |  R^-Feat  (2  ,  i)  |  +  |  R^-Feat  (  3  >  i)  |  +  |  R^-Feat ( 4 ,  i)  | 


<  a  for  any  i=l,2..k  then  go  to  2,  otherwise  proceed. 

4.  Set  k=k+l,  Feat(j,k)  =  R.  for  j=l,2,3,4. 

J 

5.  If  all  the  orientations  of  the  object  have  been  considered 
then  go  to  6,  otherwise  go  to  2. 

6.  Replace  k  by  k-1. 

In  Table  1,  a  was  chosen  to  be  0.05.  Although  there  is  no 
correlation  between  ( 14 )  and  the  error  associated  with  a  pattern 
classifier,  it  is  still  useful  for  comparisons  between  object:?. 

There  is  a  considerable  degree  of  similarity  between  the  statistics 
of  the  features  which  implies  that  views  of  one  object  can  be  derived 
from  a  linear  transformation  of  the  view  of  another  object.  On  the 
surface  this  looks  discouraging  but  may  actually  prove  to  be  of  some 
value  as  discussed  in  the  section  on  future  research. 


Initially,  the  entire  aircraft  shown  in  figure  A1  was  rotated 
through  the  significant  range  of  angles  in  order  to  calculate  the 
moments  through  using  program  three  in  the  appendix.  This 
resulted  in  a  large  computational  load  and  partial  results  indicated 
that  the  moments  were  widely  different  for  nearly  all  aspect  angles. 
Instead,  the  aircraft  was  partitioned  into  three  regions:  the  wing, 
fuselage  and  tail  as  defined  by  the  image  coordinates  in  the  appendix. 

Each  region  of  the  aircraft  was  rotated  through  the  significant 
range  of  angles  and  R  through  R^  were  calculated.  The  results  are 
summarized  in  Table  2. 

Again,  moments  R^  and  R^  are  nearly  zero  for  all  orientations 
of  the  object.  For  the  tail  section,  the  variances  of  the  moments 
are  the  smallest  implying  a  relaxed  dependency  upon  aspect  angle. 

Also,  the  moments  for  the  wing  and  fuselage  are  nearly  the  same  in 
the  case  of  R^  and  the  large  variance  on  R^  in  both  cases  leads  one 
to  conclude  that  these  two  regions  would  be  difficult  to  separate 
in  a  pattern  classifier  with  these  moments  as  features.  Evidently, 
the  tail  section  would  be  the  most  useful  for  classifying  aircraft 
type. 

Summary 

■r 

The  use  of  moments  which  are  invariant  to  a  linear  transformation 
were  investigated  for  various  objects  and  for  various  regions  of  an 
aircraft.  Examining  the  statistics  of  the  moments  it  appears  that 
it  would  be  difficult  to  classify  certain  objects  using  these  moments 
as  features  because  of  their  statistical  similarity.  On  the  other 
hand,  the  aircraft  experiments  seem  to  indicate  that  certain  regions 
of  the  aircraft  could  be  more  easily  classified  than  others. 


Suggestions  for  Future  Work 


The  original  purpose  of  this  work  was  to  minimize  the  number  of 
feature  sets  required  to  characterize  an  object,  or  a  region  of  the 
object,  with  the  intention  that  this  reduced  feature  set  could  then 
be  used  for  discriminating  between  classes  of  objects.  It  seems 
apparent  from  the  results  that  many  views  of  an  object  are  just 
linear  transformations  of  another  view  of  the  object  and  so  moments 
can  be  used  which  result  in  a  reduced  feature  set  i.e.  the  aspect 
angle  dependency  has  been  relaxed.  However,  examining  the  values 
of  R^,  R2,  and  R^  it  is  also  apparent  that  they  won't  be  able 
to  separate  classes  with  any  significant  degree  of  accuracy.  On  the 
other  hand,  there  is  an  important  conclusion  which  can  be  drawn 
from  these  data  which  is  essentially  a  summary  of  the  contribution 
of  this  investigation  subject  to  further  experimentation. 

The  computation  of  R^  through  R^  seems  to  illustrate  that  (not 
only  are  many  views  of  an  object  just  linear  transformations  of 
another  view  of  the  object  but  that)  a  view  of  one  object  can  be 
derived  from  a  linear  transformation  of  the  view  of  another  object. 
This  is  seen  from  the  similarity  between  the  statistics  of  R^  through 
R^  for  different  objects.  This  suggests  that  some  regions  are 
likely  to  cause  problems  in  a  pattern  classifier  if  features  are 
extracted  from  a  silhouette  that  includes  this  region.  Consequently, 
such  objects  or  regions  with  statistically  similar  moments(which  are 
invariant  to  a  linear  transformation)  constitute  ambiguous  shapes. 

In  the  context  of  pattern  recognition,  the  features  derived  f^om 
these  ambiguous  shapes  occur  near  the  sub-space  boundaries  and 
cause  problems  in  a  pattern  recognition  algorithm.  As  a  result, 
shape  algorithms  are  then,  by  definition,  a  collection  of  ad  hoc 
rules  needed  to  correct  erroneous  judgements  caused  by  these  so 
called  uifficult  cases.  This  concept  is  similar  to  that  taken  by 
Blesser  and  Shillman  (7,8)  in  their  theory  of  character  recognition 
in  which  ambiguously  shaped  characters  formed  the  basis  for  their 
theory.  With  this  concept  of  an  ambiguous  shape  in  mind  it  seems 
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that  further  investigations  with  region  analysis  techniques  (19) 
would  determine  whether  they  can  aid  in  identifying  these  ambiguous 
shapes.  Also,  it  has  been  established  in  a  number  of  cases  in 
character  recognition  that  techniques  such  as  polygonal  approximation 
greatly  reduce  the  influence  of  noise  on  the  classification  error  (19), 
a  problem  which  is  significant  when  using  moment  invariants. 

Actually  what  is  needed  in  this  area  is  a  theory  of  shape  based 
on  human  perception.  The  problem  with  any  shape  description  lies 
in  the  fact  that  there  exist  ambiguous  shapes.  For  example,  how 
much  shape  distortion  is  permitted  before  a  round  shape  becomes 
an  oval  shape?  Indeed,  it  is  these  ambiguous  shapes  which  show  up 
as  difficult  cases  in  a  pattern  recognition  algorithm.  Consequently, 
the  argument  presented  here  is  based  on  ambiguous  shapes  rather  than 
archetyped  shapes.  Because  objects  differ  markedly  from  the  "norm" 
shape,  a  computer  algorithm  should  measure  the  distance  between  the 
unknown  shape  and  the  boundary,  rather  than  between  the  unknown  and 
the  norm  shape.  In  effect,  the  boundary  between  classes  is  represented 
by  ambiguous  shapes.  The  problem  here,  is  to  find  that  boundary. 

To  do  so,  a  specific  knowledge  about  human  classification  of  shapes 
must  first  be  obtained. 

A  knowledge  about  human  classification  of  shapes  can  be  found 
through  psychological  experimentation.  First,  define  three  classes 
of  attributes  that  an  image  might  have  (7,8):  physical  attributes, 
perceptual  attributes  and  functional  attributes.  Physical  attributes 
are  the  parts  of  the  image  usually  described  in  geometric  or 
topological  terms  e.g.  line,  angles,  areas,  perimeters  etc.  Perceptual 
attributes  are  the  qualities  which  are  perceived  ,  for  example,  two 
lines  may  be  perceived  as  being  equal  even  though  they  are  physically 
unequal.  Functional  attributes  provide  the  description  of  the  image 
e.g.  does  a  human  subject  classify  the  subject  as  being  round  or 
oval? 

Consider  the  attribute  height  equals  width  (HEW)  which  might 
be  used  to  classify  an  image  as  round  or  oval.  The  state  of  the 
physical  attribute  is  determined  by  a  physical  measurement  which 


mathematically  answers  the  question,  does  height  equal  width?  The 
state  of  the  perceptual  attribute  is  determined  by  experimentally 
answering  the  question,  do  you  see  height  equal  to  width?  The  state 
of  the  functional  attribute  can  only  be  determined  by  answering, 
is  the  image  round  or  oval?  For  example,  a  truly  circular  image 
possesses  HEW  in  all  three  cases:  it  is  physically  present, 
perceived  as  being  present  and  also  functions  as  being  present 
since  such  an  image  would  be  labelled  as  round.  Introducing  a 
slight  distortion  in  the  circularity  of  the  image  (which  would 
cause  drastic  changes  in  the  moment  invariants)  ,  the  height  is  not 
equal  to  the  width  so  the  physical  attribute  is  not  present.  The 
perceptual  attribute  is  present  since  they  would  be  perceived  as 
being  equal  and  the  functional  attribute  would  be  present  since 
a  human  would  label  such  an  image  as  round.  Introducing  still  more 
distortion  in  the  circle  so  that  it  begins  to  approach  an  oval,  the 
physical  attribute  is  not  present,  the  height  is  perceived  as  being 
unequal  to  the  width,  but  yet  it  is  unclear  as  to  whether  or  not 
the  HEW  is  functionally  present  since  the  image  cannot  be  labelled 
either  round  or  oval  with  any  high  degree  of  confidence.  The  image 
is  therefore  an  ambiguous  image  due  to  the  transitional  state  of  the 
functional  attribute  HEW.  When  the  image  is  clearly  oval,  all  three 
attributes  are  not  present  and  the  image  would  be  labelled  as  oval. 

The  problem  now  is  to  find  the  mapping  between  the  physical  attributes 
and  the  functional  attributes  called  Physical  to  Functional 
Rules  ( PFRs ) ,  (7,8).  * 

The  PFR  for  this  example  would  involve  finding  a  and  b  such 

that 


b  <  H  <  a 

W 

H  >  a 

W 

H  <  b 

W 


Round 

oval 

oval 


Shillman  (7,8}  has  developed  several  experimental  methods  for 
finding  ambiguous  characters.  In  our  case,  the  ambiguous  image 
and  hence  the  PFR  can  be  found  through  labeling  experiments,  reaction 
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time  experiments  or  goodness  experiments. 

In  labeling  experiments,  subjects  would  label  shapes  along 
trajectories  which  go  from  one  shape  subspace  to  another.  The 
boundary  would  then  be  identified  as  the  shape  or  series  of  shapes 
which  are  assigned  the  two  shape  labels  with  equal  probability. 

In  reaction  time  experiments,  the  amount  of  time  required  to 
label  a  shape  can  be  measured.  Near  the  boundary,  subjects  spend 
a  greater  than  average  amount  of  time  identifying  certain  stimuli 
(7,8).  So  reaction  time  is  another  technique  for  determining  inter 
shape  boundary  locations. 

Goodness  experiments  are  based  on  the  concept  of  fuzzy  sets. 
Subjects  would  be  asked  to  rate  each  stimulus,  using  the  integers 
0-5,  indicating  how  good  a  representation  the  stimulus  is,  for  say, 
round  or  oval.  It  is  argued  that  these  concepts  and  experiments 
can  be  used  successfully  to  develop  a  theory  for  describing  shape. 

These  techniques  can  also  contribute  to  developing  a  theory 
for  recognizing  images  i.e.  separating  an  image  from  the  background 
Numerous  edge  detection  algorithm  exist,  but  none  are  based  on 
human  perception.  Threshold  selection  techniques  have  been  used 
as  a  basic  tool  in  image  segmentation,  but  little  work  has  been 
done  on  the  problem  of  evaluating  a  threshold  of  an  image.  Some 
authors  evaluate  thresholds  based  on  a  busyness  criterion  and  a 
descrepancy  criterion.  However,  what  proof  is  there  that  these 
two  criteria  have  any  significance  in  terms  of  human  perception? 

The  concepts  presented  here  could  easily  develop  into  a  theory 
for  thresholding  image  boundaries. 

In  a  similar  sense,  humans  have  the  ability  to  recognize  and 
classify  objects  in  the  prescence  of  noise.  Noisy  images  have 
caused  tremendous  problems  in  pattern  recognition  based  system. 
However,  the  experiments  described  earlier  could  just  as  easily 
have  been  carried  out  with  noisy  images.  Subjects,  can  be  shown 
a  trajectory  of  computer  processed  shapes  that  have  been  corrupted 
by  noise.  Labeling  experiments  would  then  automatically  aid  in 
classifying  the  noisy  shape.  This  would  then  lead  to  the  selection 
of  a  and  b  in  the  noisy  case. 

Another  major  advantage  of  the  procedure  described  here  is 
that  it  lends  itself  to  the  design  of  a  decision  tree  classifier. 
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Finding  a  and  b  amounts  to  finding  a  threshold  at  a  decision  point 
in  the  tree. 

In  this  section  it  has  been  argued  that  theories  for  image 
recognition  and  description  can  be  developed  based  on  human 
perception.  The  net  result  would  remove  the  ad  hoc  approach  that 
presently  dominates  this  area. 


Table  1 
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Ob j  ect 

Sphere 

Cube 

Cylinder 

Cone 

Ri 

2.193 

1.041 

1.109 

0.482 

"l 

-4 

1.95  10 

0.015 

0.068 

9.16  10~3 

-1.28  10-7 

0.000 

8.24  10~7 

-1.92  10-3 

°2 

3.75  10-13 

0.000 

2.23  10-10 

1.87  10-5 

R3 

-4 

-2.11  10 

0.000 

-4 

5.5  10 

-0 .0144 

..  °i  . 

2.37  10-7 

0.000 

1.04  10-6 

6.86  10'14 

1.292 

_ 

. 

1.761 

2 

a4 

8.961 

- 

- 

7.801 

k 

1 

6 

7 

38 

x5 


Table  2 


Region 

Wing 

Tail 

Fuselage 

1.14 

0.834 

1.27 

2 

°1 

0.24 

0.056 

0.35 

*2 

-4.55  10-5 

-4.95  10-5 

-8.98  10-9 

«2 

°2 

8.29  10-9 

8.47  10-9 

2.69  10-6 

*3 

-7.32  10~3 

-4.1  10~3 

-3.06  10-2 

2 

a3 

2.9 4  1 0-5 

2.03  10-5 

1.12  10"3 

r4 

15.6 

1.96 

107.00 

2 

°4 

266.0 

6.09 

1.76  10 21 
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Appendix  1 


This  appendix  contains  the  programs  and  image  coordinates  used 
to  calculate  the  moments  R  for  j=l,2..^. 

The  first  program  calculates  the  moments  for  the  sphere  des¬ 
cribed  in  a  previous  section  as  well  as  the  mean  and  variance  of 
these  moments.  With  minor  modifications  the  program  was  also  used 
to  calculate  moments  for  a  cube,  cylinder  and  cone. 

Program  2  was  used  to  rotate  various  regions  of  the  aircraft 
shapes  shown  in  figure  A1  through  the  significant  range  of  angles 
while  calculating  the  moments  from  the  solid  silhouette.  The  last 
statement  in  the  program  determines  what  data  files  are  used  and 
hence  what  region  is  under  investigation.  The  data  files  of  the 
image  coordinates  are  in  the  latter  part  of  the  appendix. 

Next  is  figure  A1  which  contains  the  ortogonal  views  of  a 
typical  aircraft.  The  boundaries  of  these  views  were  discretized 
to  form  the  data  files  of  image  coordinates. 

Program  3  was  used  to  rotate  the  entire  aircraft  and  to  calculate 
the  moments  of  the  solid  silhpuette.  This  program  includes  boundary 
tracing  algorithms  which  had  to  be  written  to  circumvent  the  complex¬ 
ities  associated  with  certain  regions  of  the  orthogonal  views  e.g. 
concavities.  Options  are  also  included  for  working  with  the  solid 
silhouette  or  just  the  boundary  of  the  silhouette  in  the  uv  plane. 

A  smoothing  algorithm  is  also  provided  to  smooth  the  boundary  before 
calculating  the  moments. 

Data  set  1  contains  the  image  coordinates  obtained  from  the 
views  of  the  aircraft  in  figure  Al.  The  x  axis  corresponds  to  the 
roll  axis  with  the  nose  in  the  positive  x  direction.  The  y  axis 
corresponds  to  the  pitch  axis  of  the  aircraft  with  the  left  wing  on 
the  positive  y  axis.  The  z  axis  is  the  yaw  axis  with  the  tail  on 
the  positive  z  axis. 


Data  sets  2,3  and  4  are  the  image  boundary  coordinates  used 
for  the  regions  of  the  airplane  designated  as  wing  (2),  fuselage  (3) 
and  tail  ( 4 ) . 
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Program  1 


DIMENSION  XY ( 1 00 , 1 00 ) , FEAT ( 4 , 1 00 ) , XMEAN ( U ) , VAR ( 4 ) , FT ( 4 , 200 ) 
C 

DELTAX=.5 

ALPHA=.05 

DEGINC=10. 

MDEG=360./DEGINC 
DEGINC=DEGINC*2. *3. 14159/360. 


HNXYZsIO. 

NPTS=40 
HNPTSsNPTS/2 
HwXYZS = i  iNXYZ  * s2 
INC=10 

STEP= DEL TAX* *2 

NXYZ=2*(HNXYZ-INC*STEP)/DELTAX+2.*INC+2 

IDIFFsHXYZ-INC 

START=-HNXYZ-STEP 

C 

DO  50  1=1,4 
XMEAN (I)=0. 

VAR(I)=0. 

50  FEAT(I, 1 )=0. 

K195=0 

K=1 

C 

PSI=0. 

THETAs-DEGINC-3 . 1 4 1 59* . 5 
CSPS=COS(PSI) 

SNPS=SIN(PSI) 

C 

DO  800  K2= 1,19 
THETA=THETA+DEGINC 

PHI=-DEGINC  *•••  •' . .  -  '  • 

CSTI 1= COS ( THET A ) 

SMTH=SIN(THETA)  ♦ 

C1=CSTH*CSPS 
C 

DO  800  K3=1, 10 
PHI=PHI+DEGIHC 
CSPH=COS(PHI) 

SNPH=SItl(PHl) 

C 

C2=SNTU*SNPH-CSTU*SHPS*CSPH 

C3=CSTII*SNPS*SNPH+SHT1I*CSPH 

C4=CSPS*CSPH 

C5=CSPS*3NPH 

C 

DO  10  14=1, NPTS 
DO  10  15=1, NPTS 
10  XY(I4,I5)=0. 

C 


non 
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X= START 

DO  600  1=1 ,NXYZ 

IF ( I . LT . INC .OR . I .GT . IDIFF )  GO  TO  601 
X=X+DELTAX 
GO  TO  602 
601  X=X+STEP 
602  Y=START 

xx=x*x 

do  6c o  j=  i,:;xyz 

IF(J.LT. INC. OR. J.GT. IDIFF)  GO  TO  603 
Y=Y+DELTAX 
GO  TO  604 
603  Y=Y+STEP 
604  Z= START 

YY=Y*Y 

DO  600  K7=1,NXYZ 

IF(K7.LT. INC. OR. K7.GT. IDIFF)  GO  TO  605 

Z=Z+DELTAX 

GO  TO  606 

605  Z=Z+STEP 

606  CONTINUE 

IF (XX+YY+Z*Z.GT.HNXYZS)  GO  TO  600 
C 

XP=C1*X+C2*Y+C3*Z+HNPTS 

YP=SNPS*X+C4*Y-C5*Z+HNPTS 

I1=XP 

I2=YP 

XP1=I1 

YP1=I2 

IF(XP-XP1 .GT.0.5)  11=11+1 
IF(YP-YP1 .GT.0.5)  12=12+1 
C 

XY(I1 , 12 ) = 1 . 

600  CONTINUE  * 

COMPUTE  MOMENTS 

u=o. 

UM10=0. 

UM01=0. 

C 

DO  700  11=1 ,NPT3 

DO  700  12= 1 ,NPTS 

IF(XY(I1 ,12) .EQ.O. )  GO  TO  700 

U=U+1. 

UM10=UM10-I1+HNPTS 
UM01=UM01+I2-(1NPTS 
700  CONTINUE 
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XBAR=UM10/U 

YBAR=UM01/U 

AA=0. 

BB=0. 

CC=0. 

A=0. 

B=0. 

c=o. 

D=0. 

DC  750  11  =  1  ,f.'PT3 
DO  750  12=1 ,NPTS 

IF(ABS(XY(I1 ,I2)).LT..1)  GO  TO  750 

TEMPX=-I 1 +HNPTS-XBAR 

TEMPY=I2-HMPTS-YBAR 

TEMP  1 =TEMPX**2 

AA=AA+TEMP1 

BB=BB+TEMPX*TEMPY 
TEMP2= ( TEMPY )  **2 
CC=CC+TEMP2 

A=A+( TEMPX) *TEMP 1 
B=B+( TEMPY )*TEMP1 
C=C+(TEMPX)*TEMP2 
D=D+( TEMPY )*TEMP2 
CONTINUE 

U2=U*U 

U3=U2*U 

R1=(AA*CC-BB*BB)/U3 

R2= ( ( A  *D-B*C ) **2) -4 .* ( A*C-B*B ) * ( B*D-C*C ) 

R2=R2/(U3*U2*U2) 

R3=AA*(B*D-CSC)-BB*(A*D-B*C)+CC*(A*C-B*B) 

R3=R3/(U3*U2) 

TEMP=2. *BB*BB-AA*CC 

ABC=AA*BB*CC 

CC2=CC*CC 

CC3=CC2*CC 

AA2=AA*AA 

AA3=AA2*AA 

R4 = A*A*CC3-6 .  *A*B*BB*CC2+6 .  *A*C*CC*TEMP 

R4  =  R4+A*D* ( 6 . *ABC-8 . *BB*BB*BB )+9 . *B*B*AA*CC2- 18 . *B*C*ABC 

R4=R4+6. *B*D*AA*TEMP+9- *C*C*AA2*CC 

R4=R4-6.*C*D*BB*AA2+D*D*AA3 

R4=R4/(U3*U2*U2) 
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C 

K195=K195+1 

FT(1,K195)=R1 

FT(2,K195)=R2 

FT(3,K195)=R3 

FT(4,K195)=R4 

C 

c 

TOT=ABS ( R 1 ) +AES ( R2 ) +ABS ( R3 ) +ABS ( R4 ) 

DO  776  K9=1,K 
C 

IF ( ( ABS( R 1 -FEAT ( 1 , K9 ) ) +ABS ( R2-FEAT ( 2 , K9 ) ) +ABS (R3-FEAT( 3 , K9 ) ) + 
1ABS(R4-FEAT(4,K9)))/TOT.LT. ALPHA)  GO  TO  790 
776  CONTINUE 

K=K+1 

FEAT(1,K)=R1 

FEAT(2,K)=R2 

FEAT(3,K)=R3 

FEAT(4,K)=R4 

790  WRITE(6 ,200)  R1 ,R2,R3,R4, K, THETA, PSI, PHI 

200  FORMAT(E10.3,2X,E10.3,2X,E10.3,2X,E10.3,I5,3F10.4) 

800  CONTINUE 
C 

C  COMPUTE  THE  MEAN  AND  VARIANCE  OF  THE  FEATURES 

C 

DO  810  1=1 ,K195 
DO  810  J=1 ,4 

810  XMEAH ( J ) rXMEAN ( J ) +FT ( J , I ) 

DO  820  J=1 , 4 

820  XMEAN ( J ) =XMEAN ( J) /K 1 95 
DO  830  I=1,K195 
DO  830  J=1 ,4 

830  VAR(J)=VAR(J)+(FT(J,I)-XMEAN(J))«2 
DO  840  1=1,4 

840  VAR  Cl) =VAR ( I ) / CK1 95— 1 ) 

WRITE(6,200)  (XMEAN(K) ,K=1 ,4) 

WRITE(6,200)  (VAR(K) ,K= 1 ,4) 

C 

STOP 

END 
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Program  2 


2^4 


DIMENSION  XXY(80,80) ,FEAT(4, 100) ,XY(200,2) ,XZ (200,2) , 
1YZ(200,2),XMEAN(4),VAR(4),FT(4,200) 

C 

SCF-.5 

READ  NUMBER  OF  BOUNDARY  POINTS 

READ (5, 35)  NBXY 
35  FORMAT ( 13) 

DO  40  1=1, NBXY 
READ (5. 45)  XY(I. 1 ) ,XY( 1,2) 

XY(I, 1 )=SCF*XY(I, 1 ) 

XY(I,2)=SCF*XY(I,2) 

40  CONTINUE 
45  F0RMAT(F5.0,F8.0) 

READ (5, 35)  NBXZ 
DO  41  1=1, NBXZ 
READ  (5,45)  XZ(I , 1 ) ,XZ (1,2) 

XZ(I, 1 )=SCF*XZ(I, 1 ) 

XZ ( I , 2 ) =SCF" XZ (1,2) 

41  CONTINUE 

READ (5, 35)  NBYZ 
DO  42  1=1, NBYZ 
READ (5, 45)  YZ(I, 1 ) ,YZ(I,2) 

YZ(I , 1 )=SCF*YZ(I, 1 ) 

YZ ( I , 2 ) =SCF*YZ (1,2) 

42  CONTINUE 


ALPHA= .05 
DEGINC=  10 .  '■ 

DEGINC=DEGINCS3. 14159/180. 
DO  50  1=1,4 
XMEAN(I)=0. 

VAR( I)=0. 

50  FEAT(I, 1)=0. 

K=  1 

K195=0 

PSI=0. 

TUETA=-DEGINC-3. 14159/2. 
CSPS=COS(PSI) 

SI1PS=SIN(PSI) 

FVW=80. 

FVWD2=-FVW/2. 

M=80 

RN=N 

MNPTS=RN/2. 

DELTAzFVW/RN 


o  o  n  non 


SEARCH  FOR  MAX  AtJD  MIN  VALUES  OF  Z 

ZMAX=FVWD2 

ZMINs-ZMAX 

XMAX=ZMAX 

XMINsZMIN 

YMAX=ZMAX 

YMINsZMIN 

DO  10  1=1  ,NBXZ 

TESTZ=XZ(I,2) 

TESTX=XZ( 1 , 1 ) 

I F  ( TESTZ . GT . ZMAX )  ZMAXsTESTZ 
IF (TESTZ . LT . ZMIN )  ZMIN=TESTZ 
IF ( TESTX . GT . XMAX )  XMAX=TESTX 
IFCTESTX.LT. XMIN)  XMIN=TESTX 
10  CONTINUE 

M= ABSC  ZMAX-ZHI:  I ) /DELTA+ 1 . 

MM= ABS ( XMAX-XMIN ) /DELTA+ 1 . 
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DO  15  1=1 ,NBXY 
TESTY=XY(I,2) 

IF (TESTY .GT .YMAX)  YMAXsTESTY 
IF ( TESTY. LT.YMIN)  YMIM=TESTY 
15  CONTINUE 

MMM=ABS ( YMAX-YMIN ) /DELTA+1 . 


WRITEC6.211)  NBXY , NBXZ , NBYZ , N 

211  F0RMATC2X, *NBXY=\I5, 'NBXZ=',I5, 'NBYZ=\I5, 'ARRAY  SIZE=',I5) 
C 

DO  800  K2= 1,19 
THETAsTHETA+DECINC 
PHI=-DEGINC 
CSTh’= COS  (THETA) 

SNTHsSIM (THETA) 

C1=CSTH*CSPS 

C 

DO  800  K3=1,10 
PHI=PHI+DEGINC 
CSPH=COS(PHI) 

SNPH=SItl(PHI) 

C 

C2=SNTH*SNPH-CSTH*SNPS*CSPH 

C3=CSTH*SNPS*SMPH+SNTH*CSPH 

C')=CSPS*CSPII 

C5=CSPS*SNPH 

C 

DO  5  1=1, N 
DO  5  J=1,N 
5  XXY(I,J)=0. 

C 

Z=ZM IN-DELTA 
DO  600  1=1, M 
Z=Z+DELTA 
Y=YMIN-DELTA 
DO  600  11=1 ,MMM 
Y=Y+DELTA 
X=XMIIJ-DELTA 
DO  600  111= 1, MM 
X=X>CELTA 


non 


CHECK  THE  SILHOUETTE  IN  THE  XY  PLANE 


L1=0 

L2=0 

L4=0 

L5=0 

DO  330  11=1 ,HBXY 

IF (ABS(XY(I1 , 1 )-X) .GT.0.001 )  GO  TO  329 
TEMP=XY( II ,2) 

IF(Y.GE.TEMP)  L1=1 
IF(Y.LE.TEMP)  L2=1 
L3=L1+L2 

329  CONTINUE 

IF(ABS(XY(I1 ,2)-Y). GT.0.001)  GO  TO  331 
TEMP=XY(I1 ,1) 

IF(X.GE.TEMP)  L4=1 
IF(X.LE.TEnr)  L5=1 
L6=L4+Lp 

331  IF(L3.EQ.2.AND.L6.EQ.2)  GO  TO  335 

330  CONTINUE 
GO  TO  600 

335  CONTINUE 

C  CHECK  THE  SILHOUETTE  IN  THE  YZ  PLANE 

L1=0 
L2=0 
L4=0 
L5=0 

DO  430  11=1 ,NBYZ 

IF(ABS(YZ(I1 , 1 )-Y) .GT.0.001 )  GO  TO  429 
TEMP=YZ(I1 ,2) 

IF(Z.GE.TEMP)  L1=1 
IF(Z.LE.TEMP)  L2=1 
L3=L1+L2 

429  CONTINUE 

IF(ADS(YZ(I1,2)-Z). GT.0.001)  GO  TO  431 
TEKP=YZ(I1 , 1 ) 

IF(Y.GE.TEMP)L4= 1 
IF(Y.LE.TEMP)L5=1 
L6=L4+L5 

431  IF(L3.EQ.2.AND.L6.EQ.2)  GO  TO  435 

430  CONTINUE 
GO  TO  600 

C  CHECK  THE  XZ  PLANE 

435  L1=0 

L2=0 
L4=0 
L5=0 

DO  440  11=1 ,NBXZ 

IF(ABS(XZ(I1 ,1)-X). GT.0.001)  GO  TO  439 
TEMP=AZU1 ,2) 

IF (Z .GE.TEMP)  L1= 1 
IF(Z.LE.TEMP)  L2= 1 
L3=L1+L2 

439  CONTINUE 

IF(AB5(XZ(I1,2)-Z). GT.0.001)  GO  TO  441 
TEHP=XZ( I 1,1) 

IF(X.GE.TEHP)  L4=1 
IF(X.LE.TEHP)  L5=1 
L6=L4+L5 

441  IF(L3.EQ.2.AND.L6.E0.2)  CO  TO  14 

440  CONTINUE 


o  n  o  o  o 


C 

GO  TO  600 
C 

c 

14  XP=-(C1*X+C2*Y+C3*Z)+HNPTS 
YP=SNPS»X+C4*Y-C5*Z+HNPTS 
I1=XP 
I2=YP 
C 

XXY(I1,I2)=1. 


600  CONTINUE 

COMPUTE  MOMENTS 

U=0. 

UM10=0. 

UH01=0. 

C 

DO  771  11=1, N 
DO  771  12=1 ,N 

IF(ABS(XXY(I1 ,12) ) .LT.0.01 )  GO  TO  771 
U=U+1. 

UM10=UM10-I1 +HNPTS 
UM0 1  =  UM0 1 +I2-HNPTS 
771  CONTINUE 
C 

XBAR=UM10/U 

YBAR=UM01/U 

C 

PA-O. 

BB=0. 

CC=0. 

A=0. 

B=0. 

C=0.  , 

D=0. 

C 

DO  750  11=1, N 
DO  750  12= 1 , N 

IF(ABS(XXY(I1, 12)). LT.0.01)  GO  TO  750 
TEMPX= - I 1 +HNPTS-XBA R 
TEMPY= I2-HNPTS-YBAR 
TEMP1=TEMPX**2 
AA=AA+TEMP 1 
C 

BB=BB+TEMPX*TEMPY 

TEMP2=TEMPY**2 

CC=CC+TEMP2 

C 

A=A+TEMPX*TEMP1 
B= B+TEMP Y  *TEMP 1 
C=C+TEHPX*TEMP2 
D=D+TEMPY*TEMP2 
750  CONTINUE 
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U2=U*U 

U3=U2*U 

C 

R1=(AA*CC-BB*BB)/U3 

R2=((A*D-B*C)**2)-4.*(A#C-B*B)*(B*D-C*C) 

R2=R2/(U3*U2*U2) 

R3=AA*(B*D-C*C)-BB*(A*D-B*C)+CC*(A*C-B*B) 

R3=R3/(U3*U2) 

TEMP=2 . *BB*BB-AA*CC 
ABC =AA  :'B5  SCC 
CC2=CC*CC 
CC3=CC2*CC 
AA2=AA*AA 
AA3=AA2*AA 
C 

R4= A*A*CC3-6 . *A*B*BB*CC2+6 . *A*C*CC*TEMP 
R4=R4+A*D* ( 6 . *ABC-8 . *BB*BB*BB)+9 . *B*B*AA*CC2- 1 8 . *B*C*ABC 
C 

R4=R4+6. *B*D*AA*TEMP+9 . *C*C*AA2*CC 

R4=R4-6.*C*D*BB*AA2+D*D*AA3 

R4=R4/(U3*U2*U2) 

C 

K195=K195+1 

FT(1,K195)=R1 

FT(2,K195)=R2 

FT(3,K195)=R3 

FT(4,K195)=R4 

TOT= AES ( R 1 ) +ABS ( R2 ) +ABS ( R3 ) +ABS ( R4 ) 

DO  776  K9=1,K 

IF ( ( ABS ( R 1 -FEAT ( 1 , K9 ) ) +ABS ( R2-FEAT ( 2 , K9 ) ) +ABS ( R3-FE AT ( 3 , K9 ) ) + 
1 AES ( R4-FEAT ( 4 , K9 ) ) ) /TOT . LT . ALPHA )  GO  TO  790 
776  CONTINUE 
C 

K=K+1 

FEAT( 1  ,K) =R1 
FEAT(2,K)=R2 
FEAT(3,K)=R3 
FEAT(4,K)=R4 
C 

790  WRITE (6, 200)  R1 ,R2,R3,R4,K, THETA, PSI , PHI 

200  FORMAT(E10.3,2X,E10.3,2X,E10.3,2X,E10.3,I5,3F10.4) 

800  CONTINUE 
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C  COMPUTE  TOE  MEAN  AND  VARIANCE  OF  FEATURES 

C 

DO  810  I=1,K195 
DO  810  J=1,4 

8 1 0  XMEAM  ( J )  =XMEAtl  ( J )  +FT  ( J ,  I ) 

DO  820  J=1 ,4 

820  XMEAM (J)=XMEAN(J)/K1 95 
DO  830  1=1 fK195 
DO  830  J=1 ,4 

830  VAR( J)=VAR( J)+(FT( J ,I)-XMEAM( J))**2 

DO  840  1=1,4 

840  VAR ( I ) =VAR ( I ) / (K1 95- 1 ) 

WRITE(6,200)  (XMEAN (K) ,K= 1,4) 

WRITE(6,200)  (VAR(K),K=1 ,4) 

C 

500  STOP 
END 

/DATA 

/INCLUDE  FUSEXY , FUSEXZ , FUSEYZ 


o  o  o  o  ooon  onooo 


Program  3 


DIMENSION  AXY( 100, 100) ,FEAT( 4, 100) ,B8(8) ,XY(400,2) ,XZ(400,2) , 
1YZ(400,2),XYS(400,2) 

IFLAG  .NE.  1  PRINTS  OUT  THE  SOLID  SILHOUETTE 

JFLAG  .NE.  1  SMOOTHS  THE  BOUNDARY 

KFLAG  .NE.  1  PRINTS  OUT  THE  SILHOUETTE  BOUNDARY 

IFLAG= 1 
JFLAG=1 
KFLAG=1 
SCF=.5 

READ  NUMBER  OF  BOUNDARY  POINTS 

READ (5, 35)  NBXY 
35  FORMAT (13) 

DO  40  1=1, NBXY 

READ (5, 45)  XY(I , 1 ) ,XY( I ,2) 

XY( I , 1 ) =SCF*XY (1,1) 

XYS(I, 1)=XY(I, 1) 

XY(I,2) =SCF*XY( 1,2) 

XYS(I,2)=XY(I,2) 

40  CONTINUE 
45  F0RMAT(F5.0,F8.0) 

READ (5, 35)  MBXZ 

DO  41  1=1 ,NBXZ 

READ  (5,45)  XZ(I, 1),XZ(I,2) 

XZ (1,1) =SCF*XZ(I, 1 ) 

XZ ( I , 2 ) =SCF*XZ (1,2) 
il  CONTINUE 

READ (5, 35)  NBYZ 

DO  42  1=1, NBYZ 

READ(5,45)  YZ(I , 1 ) , YZ (I ,2) 

YZ( I , 1 )=SCF*YZ( I, 1 ) 

YZ(I,2)=SCF*YZ(I ,2) 

42  CONTINUE 


ALPHA= .05 
DEGINC=10. 

NDEG=360./DEGINC 
DEGINC=DEGINC*3. 14159/180. 
DO  50  1=1,4 
50  FEAT( I, 1 )=0. 

K=  1 
C 

PSI=0. 

THETA=-DEGIMC-3. 14159/2. 
CSPS=COS(PSI) 

po.T  \ 


o  o  n 
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FVW=80. 

FVWD2=-FVW/2. 

N=80 
RN=N 

HNPTS=RN/2. 

DFLTA=FVW/RN 

SEARCH  FOR  MAX  AND  MIN  VALUES  OF  Z 

ZHAX=FVWD2 
ZMINs-ZMAX 
DO  10  1=1 ,MBXZ 
TESTZ=XZ(I,2) 

IF(TESTZ.GT.ZMAX)  ZMAX=TESTZ 
IFCTESTZ.LT. ZMIN)  ZMIN=TESTZ 
10  CONTINUE 

M= ABS ( ZMAX-ZMIN ) /DELTA+ 1 . 

C 

WRITEC6 ,21 1 )  NBXY , M BXZ , NBYZ , N 

211  F0RMATC2X, 'NBXY=',I5, 'NBXZ=',I5, 'NBYZ=' ,15, 'ARRAY  SIZE=',I5) 
C 

DO  800  K2= 1,19 
THETA=THETA+DEGINC 
PHI = -DEGIN C 
CSTH=COS( THETA) 

SWTHsSIMC THETA) 

C1=CSTH*CSPS 

C 

DO  800  K3=1,10 
PHI=PHI+DEGINC 
CSPII=COS(PHI) 

SNPHsSIN(PHI) 

C 

C2=SNTH*SMPH-CSTH*SHPS*CSPH 
C3=CSTH*SNPS»SNPH'+SNTH*CSPH 
C4=CSPS*CSPH  -  ' 

C5=CSPS*SNPH 

C 

DO  5  1=1, N 
DO  5  J=1,N 
5  AXY(I, J)=0. 

C 

DO  600  Isl.MBXY 
X=XY(I  1) 

IFCX.EQ. 10000.)  GO  TO  600 
JJ=I+1 

YINIT=XY( I ,2) 

IF(XY(JJ,1).EQ. 10000.)  GO  TO  615 
DO  610  J=JJ ,NBXY 
IFCX.EQ. XY(J.D)  GO  TO  620 
610  CONTINUE 


615  Y=YINIT 

YFINrYINIT 
XY( I, 1 )=10000. 

GO  TO  427 
620  YFIN=XY( J,2) 

Y=YINIT-D£LTA 

IFCYFIN.LT. YINIT)  Y=YINIT+DELTA 
XY( J, 1 )= 10000. 

425  CONTINUE 

IF(YFIN.EQ.YINIT)  GO  TO  600 
IF(YFIN.LT.YINIT)  GO  TO  426 
Y=Y+DELTA 

IF(Y.GT.YFIN)  GO  TO  600 
GO  TO  427 

426  Y=Y-DELTA 
IF(Y.LT.YFIN)  GO  TO  600 

427  Z=ZMIN-DELTA 
DO  450  K39=1,M 
Z=Z+DELTA 

C  CHECK  THE  SILHOUETTE  IN  THE  YZ  PLANE 

L1  =  0 
L2=0 

DO  430  11=1  ,NBYZ 

IFCAESC YZ(I1 , 1 )-Y) .GT.0.001 )  GO  TO  430 
TEMP=YZ(I1 ,2) 

IF(Z.GE.TEMP)  L1= 1 
IF(Z.LE.TEHP)  L2=1 
L3=L1+L2 

IFCL3.EQ.2)  GO  TO  435 
430  CONTINUE 
GO  TO  450 

C  CHECK  THE  XZ  PLANE 

435  L1=0 

L2=0 

DO  440  '11  =  1  ,NBXZ  " 

IFCABSCXZCII ,1)-X). GT.0.001)  GO  TO  440 
TEMP=XZ(I1 .2) 

IF(Z.GE.TEMP)  L 1 = 1 
IF(Z.LE.TEMP)  L2= 1 
L3=L1+L2 

IFCL3.EQ.2)  GO  TO  14 
440  CONTINUE 


C 


GO  TO  450 


o  o  n  no 


C 

14  XP=(C1*X+C2*Y+C3*Z)+HNPTS 
YP=SNPS*X+C4#Y-C5*Z+HNPTS 
1 1 =XP 
I2=YP 
C 

AXY(I1 ,12)  =  1 . 


450  CONTINUE 
GO  TO  425 
600  CONTINUE 

COMPUTE  MOMENTS 


u=o. 

UM1C=0. 

UM01=0. 

C 

DO  771  11=1, N 
DO  771  12= 1 ,N 

IF(AXY(I1 ,12) .EQ.O. )  GO  TO  771 
U=U+1. 

X=I1-HNPTS 
UM10=UM10+X 
Y=I2-HNPTS 
UM01=UM01+Y 
771  CONTINUE 
C 

XBAR=UM10/U 

YBAR=UM01/U 

C 

AA=0. 

EB=0. 

CC=0. 

A=0. 

B=0. 

C=0. 

D=0. 

C 

DO  750  11= 1 ,N 
DO  750  12=1 ,N 

IF(AXY(I1 ,12). EQ.O.)  GO  TO  750 
X=I1-HMPTS 
Y=I2-HHPTS 
TEMPX=X-XBAR 
TEMPY=Y-YBAR 
TEMP1=TEKPX*»2 
AAsAA+TEMPI 
C 

bb=bb+teupx»tempy 

TEMP2=TEMPY**2 

cc=cc+tei:p2 
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A=A+TEMPX*TEMP1 
B=B+TEMPY*TEMP 1 
C= C+TEMPX*TEMP2 
D=D+TEMPY*TEMP2 
750  CONTINUE 
C 

U2=U*U 

U3=U2*U 

c 

R1=(AA*CC-BB*BB)/U3 

R2=(  (A*D-B*C)**2)-4. *(A*C-B*B)*(B*D-C*C) 

R2=R2/(U3*U2*U2) 

R3=AA*(B*D-C#C )-EB*(A*D-B*C)+CC*(A*C-B*B) 

R3=R3/(U3*U2) 

TEMP=2. *BB*BB-AA*CC 
ABC=AA*BB*CC 
CC2=CC*CC 
AA2=AA*AA 
AA3=AA2*AA 
C 

R4=A*A*CC3-6. *A*B*BB*CC2+6. *A*C*CC*TEMP 
R4=R4+A*D*(6. *ABC-8 . *BB*BB*BB)+9 . *B*B*AA*CC2- 18 . *B*C*ABC 
C 

R4=R4+6.*B*D*AA*TEMP+9.*C*C*AA2*CC 
R4=R4-6. *C*D*BB*AA2+D*D*AA3 
R4=R4/(U3*U2*U2) 

C 

TOT=AES(R1)+ASS(R2)+ABS(R3)+ABS(R4) 

DO  776  K9=1,K 

IF((AES(R1-FEAT(1,K9))+ABS(R2-FEAT(2,K9))+ABS(R3-FEAT(3,K9))+ 
1ABS(R4-FEAT(4.K9)))/TOT.LT. ALPHA)  GO  TO  790 
776  CONTINUE 
C 

K=K+1 

FEAT(  1  ,K)=R1 
FEAT(2,K)=R2 
FEAT(3,K)=R3 
FEAT(4 ,K)=R4 
C 

DO  788  11=1 ,NBXY 
XY (11,1) =XYS( 11,1) 

788  XY(I1 ,2)=XYS(I1 ,2) 

C 

790  WRITE (6, 200)  R1 ,R2,R3,R4. K, THETA, PSI. PHI 

200  FORMAT(E10.3,2X.E10.3,2X,E10.3,2X,E10.3,I5,3F10.4) 

800  CONTINUE 
C 

tf (  tr  ar, .  fo  .  1 1  no  m  Ai 


o  o 


DO  60  1=1, N 

60  WRITEC6 ,200)  (AXY(I,J),J=1,N) 

61  CONTINUE 

C 

IF( JFLAG.EQ. 1 )  GO  TO  701 
C  SMOOTH  THE  BOUNDARY  POINTS 

NNsN-2 
C 

DO  700  1=1, NN 
DO  700  J=1 ,MN 

IF(AXY(I,J).EQ. BLANK)  GO  TO  700 

IF (AXY(I+1,J).EQ. SYMBOL)  GO  TO  700 

IF(AXY(I+1,J+1).EQ. SYMBOL)  CO  TO  700 

IF (AXY(I,J+1).EQ. SYMBOL)  GO  TO  700 

J1=0 

J2=0 

J3=0 

J4=0 

J5=0 

IF(AXY(I+2,J).EQ. SYMBOL)  J1=1 
IF ( AXY( 1+2 , J+ 1 ) . EQ .SYMBOL )  J2= 1 
IF ( AX Y( 1+2, J+2).EQ. SYMBOL)  J3=1 
IF(AXY(I+1,J+2).EQ. SYMBOL)  J4=1 
IF(AXY(I, J+2) .EQ. SYMBOL)  J5=1 
JSUM=J1+J2+J3+J4+J5 

IF( J3.EQ. 1 .AND. J4.EQ . 1 )  AXY ( 1+ 1 , J+ 1 ) =SYMBOL 
IF ( J2 . EQ . 1 . AND . J3 . EQ . 1 )  AXY( 1+ 1 , J+ 1 ) =SYMEOL 
IF( JSUM.EQ.O)  GO  TO  700 
IF( JSUM.GT. 1 )  GO  TO  700 
JUMP=J1* 1+J2*1+J3*2+J4*2+J5*3 
GO  TO  (1,2, 3), JUMP 

1  AXY ( 1+ 1 , J ) =SYMBOL 
GO  TO  700 

2  AXYCI+1 , J+1 )=SYMBOL 
GO  TO  700 

3  AX  Y(  I ,  J.+ 1 )  =SYMBOL 
700  CONTINUE 

701  CONTINUE 
C  FIND  THE  BOUNDARY 

IF(KFLAG.EQ.I)  GO  TO  401 
DO  300  1=1, H 
DO  300  J=1,N 

IF(AXY(I,J).EQ. SYMBOL)  GO  TO  310 
300  CONTINUE 


THE  FIRST  PIXEL  HAS  BEEN  FOUND  ON  THE  BOUNDARY 
310  XY( 1 , 1 )=I 
XY( 1 ,2)=J 
KSTART=2 
J1  =  1 
KSKIPrO 


C  LOOK  FORT! IE  NEXT  PIXEL  ON  THE  BOUNDARY 

295  J1=J1+1 

B8( 1 )=AXY(I , J-1) 

B8(2)=AXY(I+1,J-1) 

B8( 3)=AXY(I+1 ,J) 

B8(4)=AXY(I+1 ,  J+1 ) 

B8(5)=AXY(I,J+1) 

B8(6)=AXY( 1-1 , J+1 ) 

B8(7)=AXY(I-1 ,J) 

B8(8)=AXY(I-1, J-1) 

C 

DO  320  K=KSTART,o 
IF(K.EQ.KSKIP)  GO  TO  320 
IF (E8(K).EQ. SYMBOL)  GO  TO  325 
320  CONTINUE 

KST0P=KSTART-1 
DO  322  K= 1 ,KSTOP 
IF(K.EQ.KSKIP)  GO  TO  322 
IF(B8(K).EQ. SYMBOL)  GO  TO  325 
322  CONTINUE 

GO  TO  350 
C 

325  KSTART=K+6 
KSKIP=K+4 

IF(KSTART.GT.8)  KSTART=KSTART-8 
IF(KSKIP.GT.8)  KSKIP=KSKIP-8 
C 

GO  TO  (331, 332, 333, 334, 335, 336, 337, 338), K 


331 

XY( J1 , 1 )=I 

XY  ( J 1 ,2)=J-1 
GO  TO  329 

332 

XY ( J 1 , 1 )=I+1 
X Y( J 1 , 2 ) = J— 1 
GO  TO  329 

333 

XY  C  J 1 , 1 )=I+1 
XY ( J 1 ,2)=J 

GO  TO  329 

334 

XY(J1 ,1)=I+1 
XY(J1 ,2)=J+1 
GO  TO  329 

335 

XY( J 1 , 1 ) =1 
XY( J1 ,2)=J+1 
GO  TO  329 

336 

XY( J1 , 1 ) =1-1 
XY( J1 ,2) =J+ 1 
GO  TO  329 

337 

X Y ( J 1 , 1 )=I-1 
X Y( J 1 ,2)=J 

GO  TO  329 

338 

XY ( J 1 , 1 )=I-1 
X Y ( J 1 ,2)=J-1 

38 


C 

329  I=XY( J1 , 1 ) 

J=XY( J1 ,2) 

IF ( J 1 .GT.400)  GO  TO  400 
IF(I.NE.XY( 1,1))  GO  TO  295 
IF( J ,NE.XY( 1 ,2) )  GO  TO  235 
400  COMTIMUE 
401  CONTINUE 

DO  410  1=1, N 
DO  410  J=1,N 
410  AXY(I, J)=ELANK 
DO  420  1=1 ,J1 
11= XY (1,1 ) 

I2=XY(I,2) 

420  AXY(I1 , I2)=SYMB0L 
C 

350  WRITEC6 ,220)  J1 
220  FORMAT (2X, 15) 

C 

DO  51  1=1, N 

51  WRITE(6^2Q1 )  (AXY(I, J) , J=1 ,N) 

201  FORMAT  (2X,'128A1) 

500  STOP 
END 

/DATA 

/INCLUDE  FGR7 10XY,FGR7 10XZ , FGR7 1 OYZ , DMBFILE 
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